
#include "VertChImpl.h"
// my includes
#include "errlog.h"
#include "Mesh_grid3d.h"
#include "AlignDomain.h"
#include "grid3d.h"
#include "cxx_misc.h" // we need fixed_size_vector_with_int_codes
// standard includes
#include <vector>
#include <cmath>
#include <assert.h>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <fstream>
// Mesquite includes
#include <Vector3D.hpp>
#include <MsqError.hpp>
#include <MsqVertex.hpp>
#include <MsqMeshEntity.hpp>
#include <MeshInterface.hpp>

// using namespace std;
using Mesquite::Vector3D;
using Mesquite::MsqError;

#if defined(PRINTSTUFF) && PRINTSTUFF
#define cout_var(var)  std::cout << #var " = " << (var) << std::endl
#else
#define cout_var(var)
#endif

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*
 *                   Implementation of VertCh_base                           *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

bool VertCh_base::
choose_nodes(Mesquite::Mesh *msh, const AlignDomain *dom,
             LevelsetCode_array &ls_codes)
{
    MsqError err;

    levelset_codes = &ls_codes;
    domain = dom;

    global_patch.set_mesh(msh);
    global_patch.fill_global_patch(err);
    Mesquite::Mesh::VertexHandle *h = global_patch.get_vertex_handles_array();

    cout_var(levelset_codes->size());
    cout_var(global_patch.num_vertices());

    assert(levelset_codes->size() == global_patch.num_vertices());

#if defined(PRINTSTUFF) && PRINTSTUFF
    std::cout << "Levelset codes with node coordinates are : " << std::endl;
    for (size_t ind = 0; ind < global_patch.num_vertices(); ++ind)
    {
        Mesquite::MsqVertex &coords = global_patch.vertex_by_index(ind);
        Mesquite::Mesh::VertexHandle h = global_patch.get_vertex_handles_array()[ind];
        std::cout << std::setw(4) << ind << " : "
                << std::setw(6) << (int) (*levelset_codes)[h] << " :\t("
                << coords.x() << " , "
                << coords.y() << " , "
                << coords.z() << ") " << std::endl;
    }
#endif

    init();

    bool ret = choose();
    if (ret) // mark chosen nodes
    {
        for (size_t idx = 0; idx < PN.size(); ++idx)
            if (PN[idx].code > 0)
                (*levelset_codes)[h[PN[idx].node]] = ls_onsurf;
    }

#if defined(PRINTSTUFF) && PRINTSTUFF
    {
        size_t CHNDS = 0;
        std::cout << "Chosen nodes are :: " << std::endl;
        for (size_t idx = 0; idx < global_patch.num_vertices(); idx++)
            if ((*levelset_codes)[h[idx]] == 0)
                std::cout << std::setw(6) << idx
                    << ((++CHNDS % 10) ? " " : "\n");
        std::cout << std::endl << "Number :: " << CHNDS << std::endl;
    }
#endif

    restore();
    // check_surface_topology();


#if defined(PRINTSTUFF) && PRINTSTUFF || 1
    {
        size_t CHNDS = 0;
        std::cout << "Chosen nodes after restore are :: " << std::endl;
        for (size_t idx = 0; idx < levelset_codes->size(); idx++)
            if ((*levelset_codes)[h[idx]] == 0)
                std::cout << std::setw(6) << idx
                    << ((++CHNDS % 10) ? " " : "\n");
        std::cout << std::endl << "Number :: " << CHNDS << std::endl;
    }
#endif

    done();
    return ret;
}

/////////////////////////////////////////////////////////////////////////////
// init() allocates memory and initializes vectors

void VertCh_base::init()
{
    init_crossed_edges();
    init_pending_nodes();
}

bool operator<(const VertCh_base::CrossedEdge &L,
        const VertCh_base::CrossedEdge &R)
{
    return (L.values[0] < R.values[0]) ||
            ((L.values[0] == R.values[0]) && (L.values[1] < R.values[1]));
}

void VertCh_base::init_crossed_edges()
{
    Mesquite::MsqError err;
    CE.clear();
    vert_NCE.resize(global_patch.num_vertices(), 0);

    for (size_t el = 0; el < global_patch.num_elements(); ++el)
    {
        Mesquite::MsqMeshEntity &T = global_patch.element_by_index(el);
        // size_t nb_nodes = T.vertex_count();
        const size_t *T_nodes = T.get_vertex_index_array();
        Mesquite::EntityTopology topo = T.get_element_type();
        size_t nb_edges = Mesquite::TopologyInfo::adjacent(topo, 1);
        for (size_t i = 0; i < nb_edges; ++i)
        {
            // indices of nodes of edge local to T
            const unsigned *edge_l_nds =
                    Mesquite::TopologyInfo::edge_vertices(topo, i, err);
            // indices of nodes of edge in global_patch
            size_t edge_g_nds[2] = {
                T_nodes[edge_l_nds[0]], T_nodes[edge_l_nds[1]]
            };
            if (is_edge_crossed(edge_g_nds))
            {
                sort2(edge_g_nds);
                // try to find edge in CE
                std::vector<CrossedEdge>::iterator ce_it;
                ce_it = std::lower_bound(CE.begin(), CE.end(), &edge_g_nds[0]);
                if (ce_it == CE.end() || edge_g_nds < (*ce_it))
                {
                    CE.insert(ce_it, edge_g_nds);
                    vert_NCE[edge_g_nds[0]]++;
                    vert_NCE[edge_g_nds[1]]++;
                }
            }
        }
    }
}

void VertCh_base::init_pending_nodes()
{
    // count pending nodes
    size_t NPN = 0;
    for (size_t nd = 0; nd < global_patch.num_vertices(); nd++)
        if (is_node_pending(nd))
            ++NPN;
    // populate pending nodes;
    PN.reserve(NPN);
    for (size_t nd = 0; nd < global_patch.num_vertices(); nd++)
        if (is_node_pending(nd))
        {
            double dist;
            Mesquite::MsqError err;
            Mesquite::MsqVertex &ncoord = global_patch.vertex_by_index(nd);
            Vector3D normal, proj;
            domain->closest_point
                    (global_patch.get_vertex_handles_array()[nd],
                     ncoord, proj, normal, err);
            dist = dist_l2(ncoord.to_array(), proj.to_array());
            PN.push_back(PendingNode(nd, dist, proj));
            PN.back().NCE = vert_NCE[PN.back().node];
        }
}

void VertCh_base::restore()
{
    LevelsetCode ls;
    LevelsetCode_array::iterator Ralph;
    Mesquite::Mesh::VertexHandle *handles;
    handles = global_patch.get_vertex_handles_array();
    for (size_t idx = 0; idx < global_patch.num_vertices(); ++idx)
    {
        Ralph = levelset_codes->find(handles[idx]);
        assert(Ralph != levelset_codes->end());
        assert(handle_of(Ralph) == handles[idx]);
        if (level_of(Ralph) == ls_onsurf && try_restore(idx, ls))
            level_of(Ralph) = ls;
    }
}

bool VertCh_base::try_restore(size_t node_index,
                              LevelsetCode &level)
{
    Mesquite::MsqPrintError err(std::cout);
    std::vector<size_t> adj_nodes;
    global_patch.get_adjacent_vertex_indices(node_index, adj_nodes, err);
    Mesquite::Mesh::VertexHandle *handles;
    handles = global_patch.get_vertex_handles_array();

    /* No type here! This is return variable */ level = ls_onsurf;
    for (size_t i = 0; i < adj_nodes.size(); ++i)
    {
        size_t adj_node = adj_nodes[i];
        LevelsetCode this_ls = levelset_codes->operator [](handles[adj_node]);
        if (this_ls == ls_onsurf)
            continue;
        else if (level == ls_onsurf)
            level = this_ls;
        else
        {
            if (level != this_ls)
                return false;
        }
    }
    assert(level != ls_onsurf);
    return true;
}

void VertCh_base::done()
{
    vert_NCE.clear();
    CE.clear();
    PN.clear();
}

/////////////////////////////////////////////////////////////////////////////

void VertCh_base::dump_pending_nodes(std::ostream&S)
{
    S << "Pending nodes are : " << std::endl;
    for (size_t i = 0; i < PN.size(); i++)
        S << std::setw(4) << i << " :"
            << std::setw(4) << PN[i].node << " :"
            << std::setw(4) << PN[i].code << " :"
            << std::setw(4) << PN[i].NCE << " :"
            << std::setw(10) << PN[i].dist << " :\t("
            << PN[i].proj.x() << ", " << PN[i].proj.y() << ", " \
            << PN[i].proj.z() << ")" << std::endl;
}

void VertCh_base::dump_crossed_edges(std::ostream&S)
{
    S << "Crossed edges are : " << std::endl;
    for (size_t i = 0; i < CE.size(); i++)
        S << std::setw(4) << i << " :"
            << std::setw(4) << CE[i].values[0] << ", "
            << std::setw(4) << CE[i].values[1] << " :\t"
            << std::setw(4) << CE[i].codes[0] << ", "
            << std::setw(4) << CE[i].codes[1] << std::endl;
}

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*
 * VertCh_rec.h is a generic implementation of the recursion                 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#define _in_VERTCHIMPL_CC_
#include "VertCh_rec.h"
#undef _in_VERTCHIMPL_CC_

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*
 *                       Here we begin 3D                                    *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

class VertCh3D_ceT4 : public VertCh_rec < 3 >
{
public:

    VertCh3D_ceT4()
    {
        mSignToMark = (-1);
    }
};

class VertCh3D_T4ce : public VertCh_rec < 3 >
{
public:

    VertCh3D_T4ce()
    {
        mSignToMark = (1);
    }
};

static VertCh3D_ceT4 ceT4; // recursion by elems
VerticesChooser *VerticesChooser::DEFAULT_CHOOSER_3D = &ceT4;
//static VertCh3D_T4ce T4ce;  // recursion by edges
//VerticesChooser *VerticesChooser::DEFAULT_CHOOSER= &T4ce;

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*
 *                       Here we begin 2D                                    *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

class VertCh2D_rec : public VertCh_rec < 2 >
{
public:

    VertCh2D_rec()
    {
        mSignToMark = (1);
    }
};

static VertCh2D_rec chooser_2D;
VerticesChooser *VerticesChooser::DEFAULT_CHOOSER_2D = &chooser_2D;


